The Principle of Strain Reconstruction Tomography: 
Determination of Quench Strain Distribution from Diffraction Measurements 

Alexander M. KorsunskyQ Willem J.J. Vorster, Shu Yan Zhang, Daniele Dini, and David Latham 

Department of Engineering Science, University of Oxford 
Parks Road, Oxford 0X1 3PJ, England 

C^} Mina Golshan 

Synchrotron Radiation Source, Daresbury Laboratory Keckwick Lane, Warrington WA4 4AD, England 

(N 

O j. Jian Liu 

■ Department of Chemistry, University of Durham South Road, Durham DH1 3LE, England 

C/} \ (Dated: February 2, 2008) 

Evaluation of residual elastic strain within the bulk of engineering components or natural objects 
is a challenging task, since in general it requires mapping a six-component tensor quantity in three 
dimensions. A further challenge concerns the interpretation of finite resolution data in a way that is 
O ' commensurate and non-contradictory with respect to continuum deformation models. A practical 

solution for this problem, if it is ever to be found, must include efficient measurement interpretation 
and data reduction techniques. In the present note we describe the principle of strain tomography 
by high energy X-ray diffraction, i.e. of reconstruction of the higher dimensional distribution of 
strain within an object from reduced dimension measurements; and illustrate the application of this 
principle to a simple case of reconstruction of an axisymmetric residual strain state induced in a 
cylindrical sample by quenching. The underlying principle of the analysis method presented in this 
paper can be readily generalised to more complex situations. 
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I. INTRODUCTION 



The diffraction of neutrons and synchrotron X-rays provides a unique non-destructive probe for obtaining infor- 
mation on strains deep in the bulk of engineering components and structures. It has become a mature tool for the 
H ' determination of residual strain states in small coupons, and developments are under way to establish the facilities 
t-H ' for performing high resolution measurements directly on larger engineering components. 

There are two principal methods for extracting residual elastic strain information within objects with the help of 
diffraction, namely, the angle dispersive (monochromatic beam) technique, and the energy-dispersive (white-beam) 
technique. 

In the monochromatic mode the beam is first passed through an optical element that transmits only photons with 
the energy or wavelength lying within a narrow bandwidth. The optical element in question (the monochromator) is 
• usually a highly perfect crystal placed in the beam in reflection orientation (Bragg mode) or transmission orientation 
(Laue mode). The band pass depends on the quality of the crystal, and can be as narrow as 10~ 4 or better. In 
many practical situations this accuracy may exceed requirements, and broader band pass is in fact beneficial, since it 
increases flux on the sample, the number of crystals within a polycrystalline sample that contribute to the pattern, 
etc. 

Diffraction patterns are collected by employing a detector scanning the scattering angle 2q, or a position sensitive 
detector capable of registering total photon flux simultaneously at several positions along a line or over a two- 
dimensional surface. This mode allows accurate determination of diffraction peak intensity, shape and position. 
However, it usually requires significantly longer counting times in comparison with the white beam mode in order 
to collect the data from comparable sections of the diffraction pattern, primarily due to the reduction of flux by 
d ' monochromation, but also due to the necessity of scanning the detector. 

Energy dispersive setup allows multiple diffraction peaks to be collected simultaneously, thus achieving particularly 
efficient counting statistics at energies above 30 keV 0. The accuracy of determination of individual peak position 
and shape resolution in the white-beam mode is usually related to the resolution of the energy-dispersive detector, 
but can in fact be several orders of magnitude better £| . The accuracy of interpretation in terms of latticeparameters 



and hence strain can be significantly improved by using multiple peak analysis or whole pattern fitting 
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FIG. 1: Geometry of the incident and diffracted beams and the gauge volume 'lozenge'. 



A particular feature of high energy X-ray diffraction is that the radiation is primarily scattered forward, i.e. in 
directions close to that of the incident beam. Therefore small diffraction angles have to be used, usually 29 < 15°. 
Two difficulties immediately follow. 

Firstly, it is difficult to measure strains in directions close to that of the incident beam. This is due to the fact 
that the scattering vector is always the bisector of the incident and diffracted beams. Hence for high energy X- 
rays the strain measurement directions form a shallow cone. For a scattering angle of 29 this cone has the angle of 
(180° — 29)/ 2 = 90° — 9. In practice this means that strain directions accessible for the high energy X-ray diffraction 
techniques are close to being normal to the incident beam. This situation is in stark contrast with that encountered 
in laboratory X-ray diffraction where near backscattering geometry is used, and measured strains are in directions 
close to being parallel with the incident beam. 

Secondly, because of the small scattering angle the gauge volume becomes a strongly elongated lozenge (Fig[Q. 
The total length of the scattering gauge volume is 



where L is the length of the gauge volume, hi denotes the height of the incident beam, and hu denotes the height of 
the diffracted beam. If these two heights a both equal to h, then equation (JJJ simplifies to 



For a scattering angle 29 of 10.7° the ratio L/h is close to 10.7. As a consequence the resolution in the direction of the 
incident and scattered beam is found to be degraded about ten-fold compared to the resolution in the perpendicular 
plane. 

The motivation for the present paper lies in the desire to overcome these limitations by combining the information 
available from multiple measurements. In many respects the situation is similar to that encountered in reconstruction 
tomography, when a three dimensional structure is de-convoluted, or, rather, re-constituted from 2D projections. 
There is an important distinction between the approach proposed here and the algorithms that make use of back- 
projection, but are ignorant of the properties of the object. Here instead we propose to use a parametric model of 
the object - in the present case, the unknown strain distribution - and to carry out a minimisation procedure on the 
mismatch between the measurements and their simulation. In the present paper we deliberately choose a very simple 
distribution of residual elastic strain, since this allows us to introduce the mathematical without obfuscating it with 
too many details. However, the same fundamental approach can be applied to more complex situations, and this work 
is under way. 

The use of diffraction strain analysis for residual stress characterisation holds out the prospect of furnishing an 
important validation tool for tuning up manufacturing processes so as to reduce distortion and residual stress, and 
to increase the durability of components. In recent years powerful novel modelling approaches have been put forward 
for the analysis of residual elastic strain distributions and the corresponding residual stresses, most notably using the 
eigenstrain modelling approach Q ■ 

Heat treatment procedures play a significant role in controlling the metallurgical processes within components 
that deliver the desired alloy properties in terms of the combination of strength, hardness, toughness, ductility, etc. 
Quenching is often used in order to provide the rate of cooling required in order to form precipitation-hardened 
microstructures, to create metastable phases with particular properties, such as martensite, etc. However, invariably 
quenching leads to the creation of high thermal gradients that result in differential contraction or expansion of 
neighbouring material volumes. This, accompanied by the strong dependence of yield stress on temperature, usually 
observed in metallic alloys, leads to ready plastic deformation in the hotter parts of the specimen or component 
subjected to quenching treatment, i.e. within the bulk of material and away from the surface. As a result of the 
plastic strain distribution 'frozen' into the material, a residual elastic strain (r.e.s.) and the corresponding residual 
stress distribution arise. 



L = {hi + h D ) cos 9/ sin 29 = (hi + h D )/{2sin9), 



(1) 



L/h= l/sin9 



(2) 
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These can be beneficial or detrimental to the service properties of the component, depending on the loading 
it experiences in service. For example, the usual result of the quenching procedure is to create a state of residual 
compression near the surface, which helps to inhibit the initiation and propagation of fatigue cracks. This phenomenon 
is used in quench and induction hardening processes. 

For the present study we have chosen the residual elastic strain state in a quenched cylinder of nickel superalloy 
IN718. The residual strain state is known to be axisymmetric, and in fact is known to follow a parabola as a function 
of radial distance from the axis. These elements of knowledge help formulate a simple model for the unknown strain 
distribution, and simplify the reconstruction procedure considerably. 

II. EXPERIMENTAL 

The sample used for the experiment was made from IN718 creep resistant nickel superalloy used in the manufacture 
of aeroengine components, such as combustion casings and liners, as well as disks and blades. The composition of the 
alloy in weight percent is approximately given by 53% Ni, 19% Fe, 18% Cr, 5% Nb, and small amounts of additional 
alloying elements Ti, Mo, Co, and Al. Apart from the matrix phase, referred to as 7, the microstructure of the alloy 
may show additional precipitate phases, referred to as 7', 7", and 5. 

The primary strengthening phase, 7", has the composition Ni 3 Nb and a body-centred tetragonal structure, and 
forms semi-coherently as disc-shaped platelets within the 7 matrix. It is highly stable at 600° C, but above this 
temperature it decomposes to form the 7' N13AI phase (between 650°C and 850°C), and 5, having the same composition 
as 7" (between 750°C and 1000°C). At large volume fractions and when formed continuously along grain boundaries, 
the S is thought to be detrimental to both strength and toughness |5j. The 5 phase that forms is more stable than 
the 7" phase, and has an orthorhombic structure 

The sample was made into the form of a cylinder of a=6.5mm in diameter and was subjected to high temperature 
solution treatment for two hours at 920°C, followed by quenching axially into cold water. This treatment results in 
the creation of an axially symmetric residual elastic strain profile which varies roughly according to the parabolic 
law with the distance from the sample axis, with the surface regions of the sample experiencing compressive axial 
residual clastic strain, and its core tensile axial residual elastic strain. As follows from the above discussion of IN718 
metallurgy, we may expect to find at least three phases present in the material following this heat treatment, the 
matrix 7 phase together with 7" and S phases, and possibly some traces of 7'. 

The sample was mounted on the Euler cradle on Station 16.3 at the synchrotron radiation source at Daresbury 
Laboratory. Apart from providing the usual rotational degrees of freedom (w, x, </>) the cradle is equipped with two 
translators allowing movement in the plane of the sample support platform. 

White beam energy dispersive set up was used to collect the data. The station insertion device, wavelength shifter, 
produces a smooth white beam spectrum with useable photons having energies up to lOOkeV. Only photons with 
energies exceeding 60keV were useful in the present experiment. The diffraction spectrum was recorded using a Li- 
drifted Ge detector (Canberra) that was cooled with liquid nitrogen. A multi-channel analyser was used to bin the 
pulses from the amplifier. 

Detector calibration was carried out using the procedure described by Liu et al |4j , firstly using a radioactive source, 
then using the diffraction pattern collected from the standard NIST silicon powder sample. 

The sample was mounted as shown in Fig[21 The translation ir-axis perpendicular to the incident beam direction 
was located on the stage attached to the Euler cradle. The stage was tilted around the axis perpendicular to the 
incident beam by the angle of 5°, i.e. half the scattering angle of 10°. Slit apertures of 1.2mm for the incident beam 
and 0.1mm for the diffracted beam were used in the experiment, producing the gauge volume of the type shown in 
Figffl The symmetry axis of the cylindrical sample was designated as z-axis, and the y-axis was normal to both x 
and z directions and formed a right handed coordinate system. 

III. DIFFRACTION DATA INTERPRETATION 

In order to determine the average value of the lattice parameter within each gauge volume, Rietveld refinement is 
carried out on the diffraction patterns using GSAS. 

Our preferred data treatment procedure has been described by J. Liu et al. Q and involves careful stepwise 
calibration of the detector characteristic, followed by progressive improvement of the fit quality. Detector calibration 
begins with the use of an Am209 radioactive source that provides photon flux at three distinct energies, although 
two of these energies are relatively close to each other, approximately at 15, 20 and 60keV. This is insufficient to 
determine the non-linear detector characteristic, but allows approximate linear calibration. At the next stage a 
diffraction pattern is used that was collected from a powder of a standard material (e.g. NIST Si) with precisely 
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FIG. 2: Schematic diagram of the experimental setup. 

known lattice parameter. Starting with the previously known linear detector calibration, quadratic characteristic of 
the detector is sought by requiring the best fit to the diffraction pattern. It is important to note, in passing, that if 
the diffraction angle is not known precisely, it can be determined at this stage together with the channel-to-energy 
conversion. 

Pattern refinement began with the primary fee gamma phase in the nickel superalloy IN718. This is known to 
co-exist in coherent form with the gamma-prime phase. It was found, however, that this did not provide an adequate 
fit to the experimental diffraction pattern, particularly with the peaks at 71 and 73 keV remaining unaccounted for. 

FIgEI provides an illustration of the GSAS fit to the diffraction pattern. The use of two phases in the fit (7 and 
6) allows excellent description of the data, as shown by the difference plot (bottom curve) . The diffraction pattern 
(cross markers and continuous fit curve) is shown with the background (smooth curve) subtracted. The background 
was quite high in this experiment, since measurements were made at the higher limit of the energy range accessible 
on station 16.3 at the SRS. 

The a parameter of the majority matrix 7 phase was used in the calculation of the lattice residual elastic strain 
according to the usual formula, 



where ao denotes the reference 'strain-free' value of the lattice parameter. We do not presuppose this value, but 
instead carry the above formula in 'dynamic' form through our calculations below, i.e. retain as a parameter within 
the computation, thus being able to adjust it at the end to satisfy the conditions of equilibrium. 

IV. STRAIN RECONSTRUCTION TOMOGRAPHY 

The approach adopted in this study rests on the postulate that it is necessary to make informed assumptions about 
the nature of the object (or distribution) that is being determined in order to formulate and regularise the problem. 
Interpretation of any and every measurement result incorporates a model of the object being studied. To give a very 
general example, when a ruler is pressed against the sample to determine its length, the implication is that the sample 
and ruler surfaces are collocated all along the measured length. If that is so, then the reading from the ruler is correct; 
otherwise, however, the measurement will be a lower estimate. 

As an example of more direct relevance to the present subject, consider the situation that arises when the residual 
stress within the sample surface is 'measured' using the sin 2 ip method, a whole array of assumptions is being made, 
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FIG. 3: Illustration of the GSAS fit to the diffraction pattern. Three phases used in the fit (gamma, gamma prime and 
orthorhombic) allow excellent description of the data, as shown by the difference plot (bottom curve). The diffraction pattern 
(cross markers and continuous fit curve) is shown with the background (smooth curve) subtracted. The background was quite 
high in this experiment, since measurements were made at the higher limit of the energy range accessible on station 16.3 at the 
SRS. 

namely (i) that the material is uniform, isotropic and continuous, (ii) that strain values measured at different angles 
of tilt, tp, are obtained from the same group of grains within the same gauge volume; (iii) that the component of 
stress normal to the sample surface vanishes within the gauge volume; etc. All of the above assumptions are in fact 
approximations, or, in worst cases, entirely wrong. 

In developing the mathematical technique presented below we propose to go one step further and adopt the approach 
that it is preferable to develop measurement interpretation procedures that explicitly incorporate a model of the 
measured object. This approach would seem particularly appropriate in cases when each individual measurement 
(each data point) in fact represents a convolution of the distribution being investigated, e.g. the average value of 
residual elastic strain being measured over the diffraction gauge volume. 

In the case in hand, the unknown distribution of axial residual elastic strain (r.e.s.) is known to possess a strong 
symmetry property due to the nature of the sample geometry. Due to its axial symmetry, r.e.s. can be modelled 
using a function of radial position only. For ease of analysis and without significant loss of generality, we represent 
the unknown strain distribution as a function of radial position within the sample by the following truncated power 
series: 

n n 

e ( r ) = X! c * e »( r ) = X! CiT ^ r - a ' ( 4 ) 

»=0 i=0 

Here a denotes the radius of the cylindrical specimen, n + 1 is the total number of terms, and Cj denotes the unknown 
coefficients, yet to be determined, within the series of radial power functions, &i{r) = r\ The cylindrical polar system 
of coordinates (r, 9, z) is used, for which the z-axis coincides with the axis of symmetry of the quenched cylinder. 
For the purposes of subsequent discussion we shall also write the same expression in the Cartesian coordinate system 
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(x, y, z) which shares the z axis with the polar system, so that 



e{x,y) =^2ci(^x 2 +y 2 ) . (5) 



We now proceed to simulating the result of diffraction measurements. In the experiments described in the present 
study the incident beam height amounted to 1.2mm and the diffracted beam height to 0.1mm. From equation Q 
it follows that, for scattering angle of 28 = 10°, the length of the gauge volume therefore amounted to 7.5mm and 
spanned the full depth of the sample at all positions. It is also possible to note, from considering FigQ]in detail that 
the length of the gauge volume having uniform height is given by 

L = hi/(2siad), (6) 

and is equal to 6.9mm in the case considered. Thus, the constant height part of the scattering volume spanned the 
full depth of the section of the sample made by the plane z = const, Fig|3 

The next significant assumption made in our interpretation was that the strain value obtained from the diffraction 
profile at each point in the x-scan (Fig0J is given by the average value of strain within the sampling volume, 



JJ Ai e(r)dA 



Here residual elastic strain in the sample, e, is assumed to be a function of radial position r only, dA refers to 
the element of cross-sectional area of the sample. Area corresponds to the 'footprint' of the scattering gauge volume 
within the sample centred, in the x direction, at position Xj within the sample, 



Aj : Xj — Ax/2 < z < Xj + Ax/2, -\J a 2 - x 2 < y < \J a 2 - x 2 . (8) 

Here Ax refers to the beam width in the horizontal [x) direction that was set to 1mm in the experiment. 
Using r = y/ x 2 + y 2 , together with the definition of area Aj of equation JHJ, allows equation (7J to be used as a 
definition of the following matrix ey : 



_ _ JJ^e^x 2 +y 2 )dxdy 

JJ Ai dxdy ■ ' ] 

Due to the problem's linearity, the predicted measured strains at points Xj are given by 

n 

e(xj)=J2 c ^i- ( 10 ) 

i=0 

We now form a measure of the goodness of fit provided by the model as the functional 



K K 
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G = ^(e(x J )- e fe)) 2 =^ 5>e«-e(ay) , (11) 

3=1 3=1 \i=0 I 

where e(xj) denotes the strain measured at position Xj by diffraction. 

The search for the best choice of model can now be accomplished by minimising G with respect to the unknown 
coefficients, Cj, i.e. by solving 

grad Ci G = {dG/d Cl ) =0, i = 1, n. (12) 

Due to the positive definitcness of the quadratic form (|llfl . the system of linear equations (|12f> always has a unique 
solution that corresponds to a minimum of G. 

The partial derivative of G with respect to the coefficient Ci can be written explicitly as 

K I n \ / n K K \ 

dG/dci = 2 e.jj c^j - e(xj) = 2 X! '' k ^ r -/~^ ~'^2' s ij € ( x j) = °- ( 13 ) 

j=l \fc=0 / \k=0 3=1 3 = 1 I 
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FIG. 4: Illustration of the gauge volume positions with respect to the specimen. 

We introduce the following matrix and vector notation 

E = {%}, e={e( Xj )}, c = { Cl }. (14) 
Noting that in the context of equation 1|13JI notation ekj corresponds to the transpose of the matrix E, the entities 
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appearing in equation l|14|l can be written in matrix form as: 



K K 



A = e V e kj = EeT ' b = X! e v e ( x i) = E y- 



(15) 




Hence equation assumes the form 



A c G = 2(Ac-b) =0. 



(16) 



The solution of the inverse problem has thus been reduced to the solution of the linear system 



Ab 



c 



(17) 



for the unknown vector of coefficients c = {ci}. 

The interpretation procedure described above is most conveniently implemented in Mathematica since it allows 
expressions such as equation (J7J involving integration or averaging commands to be coded in the form of functions, 
and vectors and matrices to be generated by the application of Listable functions or using the Table command. 

We now turn to the selection of the particular strain distribution model, following the general form given in equation 
We immediately note that the linear term must always be absent from the expression in equation to avoid 
the appearance of physically unacceptable strain gradient discontinuity at r = on the cylinder axis. Due to the 
simplicity of the problem, the appropriate form of function in equation J3J describing the radial strain variation should 
be chosen in the form 



In fact, tests have been performed with larger number of terms (up to 10), but as expected they did not produce an 
improvement in the quality of approximation. 

In order to find the reference lattice spacing ag referred to in equation 0, it is most desirable to apply the axial 
(z) stress balance condition. The values of radial strain were not available from the measurements, although a twin 
detector experimental setup has been proposed to achieve this measurement simultaneously with the axial strain [j| . 
In the present case, however, quench simulations demonstrate that the radial and hoop strains are much lower than 
the axial components, amounting only to about 10% of the latter, with one component being tensile and the other 
compressive. As a consequence of this observation, the Poisson effect of transverse strains on the axial stress is not 
expected to exceed about 3%. Consequently, axial strain balance can be applied instead of stress balance, in the form 



It is possible to consider this equation as a fixed relationship between parameters Co and ci. However, we use this 
expression as an implicit equation for the unknown unstrained lattice spacing ag, thus enforcing strain balance. 



The results of the tomographic strain reconstruction are illustrated in Fig[S] by way of comparison between the 
measured diffraction strain and the modelled strain, plotted as a function of position within the x-scan centred on 
the cylindrical axis of the sample. The quality of the agreement is clearly good, if a small amount of experimental 
scatter is taken into account. In fact, it may be stated that the agreement between the model and the measurements 
is in fact the best that can be delivered by a strain function obeying a quadratic variation profile. 

FigEI shows the comparison between the reconstructed axial strain within the sample, plotted as the function of 
the radial distance from the cylindrical axis, and the modelled diffraction strain, that is plotted as a function of the 
position within the x-scan across the sample that is also centred on the cylindrical axis. The distinction between the 
two curves is immediately clear. 

Firstly, it is obvious that the axial strain is greater in magnitude than the measured diffraction strain. This 
happens because diffraction strain is in fact an average of the strain distribution within the sampling volume, and 
thus necessarily has a smaller range of variation and 'misses' the peak strain values. In the plot the axial strain 
curve is intentionally extended to the surface of the specimen, r = 3.25mm, whereas the diffraction strain curve is 
terminated at x = 2.75mm, i.e. the last point for which the gauge volume is still fully immersed in the sample. 



e(r) — cq — C2r 2 , r < a. 



(18) 




(19) 



V. 



RESULTS AND DISCUSSION 
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FIG. 5: Comparison between the measured diffraction strain (markers) and the modelled diffraction strain (continuous curve), 
as a function of position within the x-scan centred on the sample axis of symmetry. 




FIG. 6: Comparison between the reconstructed axial strain within the sample (dashed line) and the modelled diffraction strain 
(continuous solid curve), as a function of the radial position, r, and the position within the x-scan centred on the sample axis 
of symmetry. 
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The most tensile and most compressive axial strain values are both close in magnitude to about 750 microstrain. 
By way of estimate, assuming Young's modulus of about 200 GPa, the tensile stress at the axis of the cylindrical 
sample may therefore be estimated to be close to 150 MPa, and a similar magnitude compressive stress is acting at 
the surface. These stresses possess sufficient magnitude in order to affect the fatigue life of a component, and must be 
taken into account. Diffraction strain tomography provides a useful tool for their determination, and can be developed 
further to apply to more complex internal strain distributions. 
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